# 计算功率谱密度函数
import numpy as np


def glp(x, num_fft=1024):
    """
    功率谱 power spectrum
    直接平方
    :param x:
    :param num_fft:
    :return:
    """
    ps = x ** 2 / num_fft
    return ps[:num_fft // 2]




def fyp2glp(Sa, damp, w, cd):
    """
    由反应谱转换得到功率谱
    :param Sa: 加速度反应谱值
    :param damp: 阻尼比
    :param w: 圆频率
    :param cd: 场地类型
    :return:
    """
    # 根据场地类型确定外包络参数
    parameter_dict = {0: (0.3, 6, 0.5), 1: (0.6, 6, 0.45), 2: (0.8, 7, 0.35), 3: (1.2, 9, 0.25), 4: (1.6, 12, 0.15)}
    (t1, t2, c) = parameter_dict[cd]
    td_dict = {0: 7.174, 1: 7.116, 2: 8.41, 3: 10.92, 4: 15.49}
    td = td_dict[cd]  # 地震动持时取值
    M = (t2 - 0.3 * t1) / td + 1 / (2 * damp * w * td) * (
                np.exp(-2 * damp * w * t2) - np.exp(-2 * damp * w * t1)) + 1 / (c * td) * (
                    1 - np.exp(-c * (td + 0.5 * t1 - t2)))
    v = w / np.pi
    G = 4 * damp / (np.pi * w) * Sa ** 2 / (
                M * (np.sqrt(2 * np.log(v * td)) + 0.5772 / (np.sqrt(2 * np.log(v * td)))) ** 2)
    return G